package qmlt.dataset.filter;

import java.util.ArrayList;
import java.util.List;

import org.apache.commons.math.linear.Array2DRowRealMatrix;
import org.apache.commons.math.linear.RealMatrix;
import org.apache.commons.math.linear.SingularValueDecomposition;
import org.apache.commons.math.linear.SingularValueDecompositionImpl;

import qmlt.dataset.Attribute;
import qmlt.dataset.DataSetBase;
import qmlt.dataset.DataSet;
import qmlt.dataset.Instance;

public class PCAFilter implements Filter
{
    
    private int k;
    
    private double energyRate;
    
    public PCAFilter(double energeRate)
    {
        this.energyRate = energeRate;
    }

    @Override
    public DataSet filter(DataSet ds)
    {
        // dim info
        int n = ds.getInstances().size();
        int m = ds.getFeatureDefs().size();

        // prepare data
        double[][] data = new double[n][m];
        for (int i = 0; i < n; ++i)
        {
            for (int j = 0; j < m; ++j)
            {
                data[i][j] = (Float) ds.getInstances().get(i).getFeatures().get(j);
            }
        }

        // pca via svd
        RealMatrix mat = new Array2DRowRealMatrix(data);
        SingularValueDecomposition svd = new SingularValueDecompositionImpl(mat);
        RealMatrix Sigma = svd.getS();
        k = getK(Sigma);
        RealMatrix V = svd.getV();
        RealMatrix V0 = V.getSubMatrix(0, m - 1, 0, k - 1);
        RealMatrix pca = mat.multiply(V0);
        double[][] pcaData = pca.getData();

        // generate new def
        List<Attribute> newFeatureDef = new ArrayList<Attribute>();
        for (int i = 0; i < k; ++i)
        {
            newFeatureDef.add(new Attribute(Attribute.FLOAT, null));
        }

        // write data
        DataSet newDs = new DataSetBase(ds.getDef().clone(ds.getId() + "-pca"), ds);
        for (int i = 0; i < n; ++i)
        {
            Instance inst = ds.getInstances().get(i);
            Instance newInst = new Instance(inst.id, newDs);
            for (int j = 0; j < k; ++j)
            {
                newInst.getFeatures().add((float) pcaData[i][j]);
            }
            newInst.setTarget(inst.getTarget());
            newDs.getInstances().add(newInst);
        }

        return newDs;
    }

    public int getK(RealMatrix sigma)
    {
        int nrows = sigma.getRowDimension();
        int ncols = sigma.getColumnDimension();
        int dim = (nrows > ncols) ? ncols : nrows;
        
        double trace = 0;
        for (int i = 0; i < dim ; ++i)
            trace += sigma.getEntry(i, i);
        
        int k = 0;
        double s = 0;
        while (s < trace * energyRate)
        {
            s += sigma.getEntry(k, k);
            k++;
        }
        
        return k;
    }

}
